function [ normalizedMeanDiff restVar stimVar ] = collapsePowerspectrum( restPowerspectrum, stimulusPowerspectrum )
%COLLAPSEPOWERSPECTRUM Process the powerspectra of rest and stimulus sets.
%   Will take 2 wave spectra and convert them to a difference (of
%   normalized mean) spectrum and variability spectrum.

% Normalized mean
disp('Calculating means...');
meanR = nanmean(restPowerspectrum, 3);
meanS = nanmean(stimulusPowerspectrum, 3);
normalizedMeanDiff = ( (meanS - meanR)./(meanS + meanR) * 2 )';

% Just for testing stuff...
%    normalizedMeanDiff = (meanS - meanR)';
%    normalizedMeanDiff = normalizedMeanDiff / max(max(abs(normalizedMeanDiff)));

disp('Calculating variances...');
% Variance figure
restVar = nanstd(divByMean(restPowerspectrum), 0, 3)';
stimVar = nanstd(divByMean(stimulusPowerspectrum), 0, 3)';

end



function [ data ] = divByMean( powerspectrum )
%DIVBYMEAN Preprocess the data for variability

data = powerspectrum;
m = nanmean(data, 3);

for i = 1:size(data, 3)
    data(:,:,i) = data(:,:,i) ./ m;
end


end